home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_100 / 138_01 / rkst2.c < prev    next >
Text File  |  1985-03-09  |  4KB  |  179 lines
  1. /*
  2.     program:    rkst2.c
  3.     created:    03-Nov-83
  4.     by:        A. Skjellum
  5.  
  6.     modified:    14-Nov-83
  7.     by:        M. J. Roberts
  8.     and:        5-Dec-83
  9.     by:        M. J. Roberts
  10.  
  11.     modified:    25-Jul-84
  12.     by:        A. Skjellum
  13.  
  14.     Copyright 1983, 1984 (c) California Institute of Technology.
  15.     All rights Reserved.  This program may be freely distributed
  16.     for all non-commercial purposes but may not be sold.
  17.  
  18.     purpose:    illustrate the use of RKS program
  19.     update:        to test the rk4n() subroutine using a system
  20.             of two differential equations.
  21.  
  22.     uses:        rk4n() (rks.c)
  23.     summary:
  24.  
  25.         integrates the differential equation system:
  26.  
  27.         u1'(t) = 8u2(t)            u1(0) = 10
  28.         u2'(t) = 2u1(t)            u2(0) = 7
  29.  
  30.         for which the exact solution is known to be:
  31.  
  32.         u1(t) = 12exp(4t) - 2exp(-4t)
  33.         u2(t) =  6exp(4t) +  exp(-4t)
  34.  
  35.     
  36. */
  37.  
  38. /* constants */
  39.  
  40. #define SYSIZE    2    /* number of functions in system */
  41. #define    Y1ZERO    10.0    /* initial value for first equation */
  42. #define Y2ZERO    7.0    /* initial value for other equation  */
  43. #define    TSTART    0.0    /* starting time for integration */
  44. #define    TEND    1.0    /* ending time for integration */
  45. #define    STEPS    50    /* 50 steps in integration */
  46.  
  47. /* variables external to all functions */
  48.  
  49. double wvalue[STEPS+1][SYSIZE];
  50. double yarray[STEPS+1][SYSIZE],tarray[STEPS];
  51.     /* integrated solution stored here */
  52.  
  53.  
  54. /* subroutines: */
  55.  
  56. /* exact(): returns exact solution value, given t */
  57.  
  58. double exact(n,t)
  59. int n;        /*  which equation is it?  0 or 1?  */
  60. double t;
  61.  
  62. {
  63.     extern double exp();    /* exponential function */
  64.  
  65.     /*  This must find solutions for both U1 and U2.  The
  66.         exact solutions are given in the header comments to
  67.         this program, above.  */
  68.  
  69.     switch (n)
  70.     {
  71.         case 0:
  72.             return(12*exp(4*t) - 2*exp(-4*t));
  73.         case 1:
  74.             return( 6*exp(4*t) +   exp(-4*t));
  75. }    }
  76.  
  77.  
  78. /* fn(j,i,t,y):  return f(t,y) given t,y values */
  79.  
  80. double fn(j,i,t,y)
  81. int j;
  82. int i;
  83. double t;
  84. double (*y)();
  85. {
  86.     switch (i)
  87.     {
  88.         case 0:
  89.             /*  u1'(t) = 8 * u2(t)  */
  90.             return(8*(*y)(1,j,i));
  91.  
  92.         case 1:
  93.             /*  u2'(t) = 2 * u1(t) */
  94.             return(2*(*y)(0,j,i));
  95.     }
  96. }
  97.  
  98.  
  99.  
  100. /* store():  the routine to store away the W values for later reference */
  101.  
  102. double store(row, col, value)
  103. int row, col;        /*  location to store the value  */
  104. double value;        /*  the actual value to store  */
  105. {
  106.     wvalue[row][col]=value;
  107.     return (value);
  108. }
  109.  
  110. /* source():  return the  W value referenced by input parameters  */
  111.  
  112. double source(row,col)
  113. int row, col;        /*  location to look up */
  114.  
  115. {
  116.     return (wvalue[row][col]);
  117. }
  118.  
  119. /* solutn(): print solution step at console */
  120.  
  121.  
  122. solutn(j,i)
  123. int j,i;    /* element numbers  */
  124. {
  125.     printf("\nt=%7.3e, y%1d=%7.3e, y%1d_exact=%7.3e, diff=%7.3e",
  126.         tarray[j],i,source(j,i),i,exact(i,tarray[j]),
  127.         source(j,i) - exact(i,tarray[j]));
  128. }
  129.  
  130. /* main program: */
  131.  
  132. main()
  133. {
  134.     /* external declarations */
  135.  
  136.     double store(), source();
  137.  
  138.     double fn();    /* ensure that this is typed as double */
  139.  
  140.     /* local variables: */
  141.  
  142.     register int i,j;
  143.  
  144.  
  145.     double init[SYSIZE];        /* initial condition matrix */
  146.  
  147.     /* begin code: */
  148.  
  149.      printf("\n\nrkst2.c    as of 25-Jul-84\n\n");
  150.     printf("    Integrates the differential equation system:\n\n");
  151.     printf("    u1'(t) = 8u2(t)            u1(0) = 10\n");
  152.     printf("    u2'(t) = 2u1(t)            u2(0) = 7\n");
  153.     printf("    for which the exact solution is known to be:\n\n");
  154.     printf("    u1(t) = 12exp(4t) - 2exp(-4t)\n");
  155.     printf("    u2(t) =  6exp(4t) +  exp(-4t)\n\n");
  156.     printf("        for     t = %7.3e to %7.3e, with %u steps\n\n",
  157.         TSTART,TEND,STEPS);
  158.  
  159.     /*
  160.         integrate the answer from t = 0 to t = 10 sec 
  161.         STEPS points.
  162.     */
  163.  
  164.     init[0] = Y1ZERO;    /* set up initial condition matrix */
  165.     init[1] = Y2ZERO;
  166.  
  167.     rk4n(fn,source,store,SYSIZE,TSTART,TEND,STEPS,init,tarray,yarray);
  168.             /* compute the answers for 1 function */
  169.  
  170.     /*  Print out solutions  */
  171.  
  172.     for(i=0;i<SYSIZE;i++)
  173.         for(j=0;j<STEPS;j++)    /* print solution */
  174.             solutn(j,i);
  175.  
  176.  
  177.     printf("\n\nEnd of execution\n\n");
  178. }
  179.